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Low-Frequency,  Bottom-Interacting  Pulse 
Propagation  in  Range-Dependent  Oceans 

MICHAEL  D.  COLLINS 
(Invited  Paper) 


Abstract — An  asymptolk-numcrkal  model  for  low-frequency,  bot¬ 
tom-interacting  pnise  propagation  in  the  ocean  U  derived.  The  model 
works  in  the  time  domain  using  an  approach  analogous  to  the  parabolic 
equation  method  that  is  commonly  used  in  the  frequency  domain.  The 
model  handles  depth  and  range  variations  in  tbe  speed  of  sound,  den¬ 
sity,  and  attennation.  The  attenuation  is  assumed  to  depend  linearly  on 
frequency  in  the  sediment.  The  accuracy  of  the  model  is  demonstrated 
with  a  benchmark. 

Hfeywonfs— parabolic,  pulse,  attennation,  density. 

1.  Introduction 

RANGE-DEPENDENT  propagation  problems  can  be 
solved  efTiciently  with  the  parabolic  equation  (PE) 
method.  It  is  most  useful  for  low-frequency  problems  be¬ 
cause  the  number  of  operations  required  for  the  numerical 
solution  of  the  PE  increases  with  the  frequency  squared  for 
two-dimensional  problems.  The  PE  method  was  first  applied 
to  underwater  acoustics  by  Tappert  [1],  who  also  provided  a 
historical  account  of  the  PE  method.  For  this  application,  the 
PE  method  has  been  highly  developed.  In  its  present  state,  it 
can  handle  density  variations  [2]-[4],  range  dependence  [5], 
three-dimensional  variations  [6]-[8],  and  wide-angle  propa¬ 
gation  [9],  [10].  The  accuracy  of  the  PE  method  has  been 
demonstrated  with  many  benchmark  comparisons  using  nor¬ 
mal  mode  and  other  results  [11]. 

Pulse  propagation  problems  can  be  solved  with  frequency- 
domain  methods  by  solving  for  one  frequency  at  a  time  and 
summing  the  results.  The  PE  method  [12],  [13],  the  fast-field 
program  [14],  and  the  WKB  method  [15],  [16]  have  been  used 
in  this  fashion.  These  approaches  can  probably  be  improved 
substantially  for  broadband  pulses  by  working  in  the  time  do¬ 
main  because  a  significant  amount  of  effort  must  be  devoted  to 
managing  the  component  frequencies  and  performing  the  sums 
with  the  frequency-domain  approach.  Since  range-dependent 
problems  are  frequently  of  interest,  a  time-domain  approach 
analogous  to  the  PE  method  would  be  very  useful.  Two  such 
models  exist,  the  nonlinear  progressive  wave  equation  (NPE) 
[17],  which  is  an  initial  value  problem  in  time,  and  the  time- 
domain  parabolic  equation  TDPE  [l8]-[20],  which  is  an  initial 
value  problem  in  range.  In  contrast  to  the  PE  method,  these 
models  have  not  been  fully  developed. 
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was  supported  by  ONR  and  NORDA. 
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The  need  for  the  development  of  bottom  interaction  capabil¬ 
ity  into  these  modes  was  the  motivation  for  the  present  study. 
A  robust  numerical  solution  was  developed  for  the  NPE  [17], 
and  it  seems  more  natural  to  march  a  solution  of  the  wave 
equation  in  time  rather  than  range.  Thus,  we  work  with  the 
NPE  rather  than  the  TDPE.  Since  the  nonlinear  term  in  the 
NPE  is  dropped,  we  do  not  adhere  to  the  terminology  for 
this  model  and  refer  to  the  linear  version  of  the  NPE  as  the 
progressive  wave  equation  (PWE).  The  number  of  operations 
required  for  the  numerical  solution  of  the  PWE  increases  with 
the  frequency  cubed.  Thus,  this  model  is  best  suited  to  handle 
low-frequency  problems. 


II.  Derivation  of  the  PWE 


We  neglect  attenuation  for  now  and  begin  with  the  following 
wave  equation  due  to  Bergmann  [2 1  ]  for  the  acoustic  pressure 


P : 


P  dt^ 


(1) 


where  t  is  time,  c  is  sound  speed,  and  p  is  density.  Azimuthal 
dependence  is  ignored,  and  cylindrical  coordinates  are  used 
with  z  being  the  depth  below  the  ocean  surface  and  r  being 
the  horizontal  distance  from  a  point  sound  source  at  z  =  Zo- 
The  dimensions  of  the  independent  variables  have  been  re¬ 
moved  using  the  factors  Cq  =  c(zo),Po  =  p(Zo),  and  the 
characteristic  time  scale  r  of  the  source  function  F(t).  Cylin¬ 
drical  spreading  applies  for  large  r  because  the  ocean  acts 
as  a  waveguide.  Thus,  spreading  can  be  handled  by  defining 
p  =  r^'^P.  We  let  ro  ►  1,  assume  |  bp/br  |  <  |  bp/bz  |,  and 
approximate  ( 1 )  for  r  >  Tq  by 


b^p  _  1  ^  ^ 

bz^  P  bz  bz  br^  bF. 


Fallowing  the  approach  used  to  derive  the  TDPE  and  the  NPE, 
we  introduce  a  reference  frame  moving  outward  at  the  refer¬ 
ence  speed  Co  with  the  new  independent  variable  s  =  r  -  Cgt 
and  new  dependent  variable 


u(s,z,t)  =  Pis  -I-  Cot,Zj). 


(3) 


With  these  definitions,  (2)  becomes 


b^u  I  bp  bu  /  cl\b^u 

bz^  p  dz  dz  ^  \  c  V  ds- 

2co  b^u  _  1  b'u 
^  bsbt  c’  dt^ 
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The  PE  method  is  based  on  the  small-angle  asymptotic 
limit.  Long-range  propagation  in  the  ocean  results  in  nearly 
vertical  wavefronts  due  to  the  fact  that  wide-angle  rays  arc 
not  trapped  by  the  oceanic  waveguide.  In  other  words,  energy 
propagation  is  limited  to  within  the  angle  $  of  the  horizontal, 
where  t  =  tan^  <  1 .  For  the  plane  wave 

u(s,zj)  =  /lscos<<>  -1-  z  sin  (/»  -I-  (cos</>  -  l)cot]  (5) 


The  coefficient  function  U,  which  depends  weakly  on  r  and  z, 
is  ignored.  The  loss  operator  is  valid  to  leading  order  because 
it  is  correct  in  the  limit  $  0. 

In  the  time  domain,  the  Fourier  transform  is  applied  to 
decompose  the  signal  as  follows: 

u(s,z,t)  =  f  u(k,z,t)eKp(iks)dk  (11) 


with  1^1  <  4>,  we  observe  that  du/dt  =  0(t)  and  du/dz  ~ 

We  assume  that  c  =  Cq -I- 0(e),  p(r,  z)  =  p'(r^,ze'^^),  r  > 
To  =  0(€~'),  and  c^r/d  =  where  p'  is  inde¬ 

pendent  of  6,^  <  (,  and  d  is  the  ocean  depth.  Using  the 
scales  suggested  by  the  above  discussion,  we  assume  that 
u{s,  z,  0  =  u'(s,  te),  where  u'  is  independent  of  e.  TTiis 
results  in  du/dt  =  0(e)  and  du/dz  =  0(r''^^).  An  analogous 
scaling  is  commonly  used  in  the  frequency  domain  for  deriv¬ 
ing  the  PE  [1],  [4].  To  leading  order  in  €,  the  approximation 
made  i.T  going  from  (1)  to  (2)  is  valid  and  (4)  becomes 


^ 

2  \dz^ 


1  dp  du\ 
p  dz  dz) 


+  (c  -  Co) 


d^u  d^u 
ds^  ^  dsdt 


=  0. 


(6) 


We  assume  u  has  compact  support  in  s  for  all  time  and  inte¬ 
grate  (6)  with  respect  to  s  to  obtain  the  PWE 


du 

Tt 


^  r 

2  J,  \az^ 


1  dp  du\ 
p  dz  dz  ) 


ds'  -f-  (Co  -  c) 


du 
ds  ' 


(7) 


1  f  °° 

u(k,z,t)  =  ^  J  u{s,z,t)cxp{-iks)ds.  (12) 

It  follows  from  (1?)  that  3m/3/  =  0(e)  and  du/dz  = 

In  (II),  the  field  is  written  as  a  superposition  of  plane  waves 
each  multiplied  by  a  function  depending  weakly  on  t  and  z- 
A  discrete  loss  operator  analogous  to  (10)  would  act  over  the 
time  increment  At  as  follows; 

u(s,  z,t)  \  u(k,  z,  0  exp(iks  -  rj0  j  k  j  coAt'idk 

J_oe 


as  t  -*  t  +  At.  (13) 


Once  again,  the  coefficient  function  of  the  plane  wave  is  ig¬ 
nored,  and  the  loss  operator  is  valid  to  leading  order  because 
it  is  correct  in  the  limit  $  0.  Substituting  (12)  into  (13) 

and  interchanging  the  order  of  integration,  we  obtain 


u(s,  Z,  t) 


■n&CpAt  p  U{s',z,t)  ds' 

T  J_*  (T)j3CoAO^  +  (5'  -s)^ 


III.  The  Loss  Term 


as  t  t  +  At.  (14) 


Let  7/  =  (40irlog,oe)"'  and  0  be  the  attenuation  in  decibels 
per  wavelength  with  =  0(().  For  the  circular  frequency  w, 
the  complex  wavenumber  K  -  k(l  +  ii\0)  is  assumed  in  the 
sediment,  where  the  0(1)  wavenumber  k  =  oi/c.  This  for¬ 
mulation  is  used  to  model  attenuation  that  increases  linearly 
with  frequency,  which  is  in  agreement  with  experimental  re¬ 
sults  [22]-[24]  involving  various  materials  and  frequencies. 
The  linear  dependence  has  been  challenged  [25],  [26].  How¬ 
ever,  the  model  we  derive  can  be  modified  to  handle  a  general 
frequency  dependence. 

The  dependent  variable  p  in  (2)  is  replaced  with  U,  which 
is  defined  by 

p{r,z,t)  =  Uir,z)expiikoS)  (8) 

with  the  assumption  U(r,z)  —  U'(rf,Zi''^),  where  ko  = 
w/C()  and  U'  is  independent  of  e.  Substituting  (8)  into  (2)  and 
retaining  leading-order  terms,  we  obtain  the  following  PE, 
which  is  valid  for  r  >  Tq  : 

dr  2ko  \  dz^  p  dz  dz  ) 

+  i{k  -  ko)U -r,0koU.  (9) 

In  the  absence  of  the  other  operators  on  the  right-hand  side 
of  (9),  the  loss  operator  -Ti0kQU  acts  over  the  range  incre¬ 
ment  Ar  as  follows; 


A  continuous  loss  operator  L  is  obtained  by  taking  the  limit 
in  (14)  to  obtain 


A/->0l  T  J_o 


u{s’,z,t)  ds' 


Lu  =  lim  - 

A/->0  T 


i: 


(t/iScoAO^  +  {s'  -  s)^ 

u{s,  z,  Ql 

A/  J 

u{s',z,t)  -  «(s,z,0 


(t;/3CoAf)2  +  {s'  -  s)2 
In  going  from  (15)  to  (16),  we  have  used  the  identity 


(15) 
ds'.  (16) 


ij/3coAt 


ds' 


(ijjScoAO^  +  (s'  -  5)2 


=  1. 


(17) 


The  limit  in  (16)  exists  as  the  Cauchy  principal  value  of  the 
integral  with  At  -  0. 

Adding  Lu  to  (7),  we  obtain  the  PWE 


du 

Yt 


1  dp  du\ 
p  dz  dz  / 


ds'  +  (Co  - 


u{s',zj)  -  u{s,z,t)  .  , 
{s'  -  5)2 


du 

^^Ts 


(18) 


U{r.z)exp{ikoS)  -•  U{r,z)exp(ikoS  -  ri0koAr) 

as  r  r  +  Ar.  (10) 


This  PWE  differs  from  the  NPE  of  [17]  by  the  presence  of  the 
density  and  attenuation  terms  and  the  absence  of  the  nonlinear 
and  spreading  terms.  The  integral  operator  appearing  in  (18) 
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is  of  the  Kramers-Kronig  type  [27],  The  approach  we  have 
described  can  be  generalized  to  handle  an  arbitrary  complex 
dispersion  relation  by  replacing  the  function  {  A:  |  in  (13) 
by  a  complex  function  of  k. 


Substituting  (24)  and  (25)  into  (26)  we  obtain 

_ 

®  ^  =  0  (27) 


rv.  Numerical  Solution  of  the  PWE 

The  domain  is  discretized  with  grid  spacings  As,Az,  and 
Af.  The  ocean  surface  at  z  =  0  is  assumed  to  be  pressure 
release,  and  the  signal  is  assumed  to  be  trapped  within  P  = 
(0.  ^m..)  X  (0,  z„..)  due  to  the  moving  frame  and  the  lossy 
bottom.  Thus,  the  boundary  condition  u  =  0  is  imposed  on 
dP.  The  ocean  depth  is  assumed  constant  over  P.  This  is  a 
valid  leading-order  approximation  if  range  dependence  is  suf¬ 
ficiently  gradual. 

We  solve  (18)  numerically  using  the  method  of  alternating 
directions,  which  was  used  to  solve  the  NPE  [17]  and  requires 
numerical  methods  for  each  of  the  following: 

du  du 

-+(c-Co)^=0  (19) 


du  _  Co  f  “  /  ^  w  ' 

dt  2  \dz^  p  dz  dz )  ^ 


(20) 


du  v/3co  f  u(s’,z,t) -u(s,z,t)  .. 

^  =  -  1  - Ti - 7^ -  as  .  (21) 

dt  r  J_oo  (5-5)2 

The  Lax-Wendroff  scheme  [28]  may  be  applied  to  solve  the 
first-order  hyperbolic  equation  (19).  The  Courant-Friedrichs- 
Lewy  (CFL)  condition, 

max  I  c  -  Co  I  A/  s  A5  (22) 


must  be  satisfied  for  stability. 

With  the  imposed  boundary  conditions,  (20)  is  equivalent 
to 


d^u  I  dp  du  2  d^u 
dz^  p  dz  dz  ^  Co  dsdt 


(23) 


Galerkin’s  method  with  linear  test  functions  is  applied  to  dis¬ 
cretize  depth  dependence.  The  depth  grid  points  are  defined  by 
z,  =  /Az.  The  basis  functions  ♦/(z)  vanish  for  |  z  -  z,  1  >  Az, 
increase  linearly  from  0  to  1  over  z,-i  <  z  <  z,,  and  de¬ 
crease  from  1  to  0  over  z,  <  z  <  Zi+i- 
The  basis  functions  can  be  used  to  approximate  a  function 
by  a  piecewise  linear  function  with  exact  agreement  at  the 
grid  points.  We  define  m,(5,0  =  u(s,Zi,t)  and  p,  =  p(z,) 
and  obtain 


M(5,z,0sS  Ui(s,t)'*i(z)  (24) 

i 


Piz)  =  Yi  Pi^/(z)-  (25) 


Galerkin’s  method  discretizes  depth  dependence  in  (23)  by 


using  (24)  and  (25)  and  requiring  that  the  following  hold  for 
all  i : 


I  dp  du  2  a2„i 
-  ^  dz  =  0 

p  oz  oz  Co  dsdt 


(26) 


where  the  vector  u(s,  t)  contains  the  values  of  u  at  the  depth 
grid  points  and  A  and  B  are  tridiagonal  matrices  with  entries 

=  2  —  p,_i<7,_i/2  (28) 


^I.I  —  P/-l<^;-I/2  +  P/+I<3/+1/2  ~4  (29) 


■^/./4-i  —  2  —  p/+ia/+i/2 

(30) 

f'  f 

—  if  Pi  =  P,4-l 

Pi 

1  log(p,^,)-log(p,)  ifp, 

Pi  4- 1  “  Pi 

(31) 

"  3co 

(32) 

«  (2Az)' 

"  3co 

(33) 

R  - 

(34) 

We  define  the  grid  points  Sj  =  y'Ar  and  t„  =  nAl  and 
»j.n  =  u(Sj,tn)-  Integrating  (27)  over  (sj,Sj+i)  x  (f„,/„  +  ,) 
using  the  trapezoid  rule  where  necessary,  we  obtain 


AsAt 

4 


(aj  4- 1  ,n  4- 1  4- 1  .n 

+  B(Uj+i,„+i  -  Uj 


-f  Uj  „) 

4-i,/i  ~  -f  Uj  „)  =  0.  (35) 


In  the  case  c  ■  Co,  p  ®  po,  and  /3  =  0,  (18)  reduces  to  (20). 
In  this  situation,  energy  flows  from  s  =  5mM  to  5  =  0  due  to 
geometric  dispersion.  This  physical  consideration  implies  that 
(35)  must  be  applied  by  sweeping  from  5  =  5nuix  to  5  =  0, 
which  results  in  the  scheme 

^ajn+l  —  bjn^\  (36) 


M  =  B- 


AsAt  ^ 

~r-A 


(37) 


+  B(M;4-I.«4.1  -  «y+l.«  +  Uj.„).  (38) 

The  loss  operator  is  defined  by  (14)  as  At  -►  0.  Thus, 
(14)  is  used  to  solve  (21).  The  CFL  condition  (22)  and  ac¬ 
curacy  considerations  result  in  5  =  ri0CoAt/As  <  l.fln  the 
example  below,  6  =  0.023.)  Thus,  the  kernel  in  (14),  which 
approximates  a  delta  function,  is  nearly  singular  at  5'  =  5. 
TTie  integral  over  (5,  -  As/2,  5,  -(-  As/2)  is  approximated 
by  replacing  m(5,  z,  t)  with  m(5„  z,  t)  and  integrating  the  kernel 
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analytically  to  obtain 


UiSj.Z.t) 


u(s,,z,0  tan  - - — 2^^ 


Fourier  decomposition  for  approximating  broadband  source 
functions  because  the  Gaussian  itself  is  a  broadband  source 
function. 


VI.  PWE  Benchmark 


For  1  j,  the  terms  in  the  sum  are  of  order  b/(i  -  jf. 
Thus,  the  sum  may  be  truncated  making  the  loss  operator 
numerically  efficient. 


V.  The  Initial  Field 


The  PWE  is  not  valid  near  the  sound  source  because  it  is 
derived  from  a  singular  perturbation  problem.  Both  du/dt  and 
1/r  fail  to  be  0(e)  in  a  boundary  layer  near  /■  =  0.  In  the 
boundary  layer  and  within  the  angle  4'  of  the  horizontal,  rays 
are  not  significantly  affected  by  the  weak  refraction  of  the 
ocean.  Thus,  the  PWE  may  be  initialized  near  the  source  at 
t  =  to  with 


P{r,z,t) 


where  d\.  -  +  (z  ±  Zo)^  ■ 

This  solution  is  valid  for  small  r  and  e  because  refraction  in 
the  ocean  is  negligible  over  short  ranges,  and  the  propagating 
rays,  which  travel  within  the  angle  ^  from  horizontal,  do  not 
intersect  the  ocean  bottom  near  the  source.  It  breaks  down 
for  large  r  because  refraction  in  the  ocean  is  significant  over 
large  ranges,  and  propagating  rays  from  the  source  eventually 
reflect  from  the  ocean  bottom. 

The  numerical  method  used  to  solve  the  PWE  requires  a 
continuously  differentiable  source  function  with  compact  sup¬ 
port.  A  source  function  that  approximates  a  delta  function 
would  be  useful  because  it  can  be  used  to  approximate  the 
impulse  response  of  the  ocean.  The  Gaussian  source  function 
G,(t)  =  exp  [-(vt)^],  which  was  used  in  [17],  has  these 
properties. 

Since  the  Gaussian  source  function  approximates  a  delta 
function,  it  is  useful  for  convolution.  A  given  source  function 
may  be  approximated  by 


FU)  =  '2  -t")  (41) 

n 

provided  v  is  sufficiently  large.  The  coefficients  a"  are  de¬ 
termined  by  requiring  the  approximation  to  be  exact  at  the 
points  which  do  not  necessarily  correspond  to  the  grid 
points  t„  and  may  be  spaced  irregularly.  This  constraint  gives 
the  following  system  of  equations; 


m-¥q 

f(t"')  =  S  a"G,(r  - 1").  (42) 

The  matrix  corresponding  to  this  system  of  equations  is  re¬ 
duced  to  the  29  -f-  1  diagonals  centered  about  the  main  diag¬ 
onal,  where  9  is  a  small  integer,  due  to  the  fact  that  G,(t) 
decays  rapidly  away  from  t  =  0.  In  some  situations,  the 
approximation  given  by  (41)  might  be  more  useful  than  the 


The  Gaussian  source  C,(0  with  v  =  150  s' '  is  placed 
at  Zo  =  75  m  in  an  ocean  in  which  c  =  Cq  =  1500  m/s. 
Ocean  depth  is  200  m  for  r  <  4  km,  linearly  decreasing  from 
200  to  50  m  over  4  km  <  r  <  8  km,  and  50  m  for  r  >  8 
km.  In  the  sediment,  c  =  1600  m/s,  p  =  1.5  g/cm^,  and 
|8  =  0.5  dB/X.  The  grid  spacings  are  At  =  2.5  m/co.  As  =  1 
m,  and  Az  =  1.5  m  wiA  =  400  m  and  Zmai  =  300 
m.  An  absorbing  layer  150  m  thick  is  added  below  z  =  Zmax 
over  which  j3  increases  to  prevent  reflections.  The  grid  f  is 
initialized  with  the  leading  Gaussian  front  at  r  =  350  m  and 
s  =  300  m. 

Fig.  1  contains  a  sequence  of  contour  plots  of  the  acoustic 
pressure.  The  ocean  bottom  is  marked  with  a  solid  horizontal 
line.  Solid  contours  represent  p  >  C;  dashed  contours  rep¬ 
resent  p  <  0.  The  leading  front  is  located  about  100  m  to 
the  left  of  the  leading  end  of  the  grid.  The  reflection  coeffi¬ 
cient  at  the  ocean  surface  is  - 1 .  The  reflection  coefficient  at 
the  ocean  bottom  is  approximately  - 1  for  small-angle  prop¬ 
agation.  Thus,  the  fronts  that  trail  the  leading  ffont  and  have 
reflected  from  the  ocean  surface  and  bottom  are  composed 
of  alternating  solid  and  dashed  contours.  The  signal  is  essen¬ 
tially  confined  to  F,  which  was  assumed  in  the  derivation  of 
the  numerical  solution. 

The  leading  front  is  accompanied  below  the  interface  by  an 
evanescent  wave.  This  is  analogous  to  the  evanescent  portions 
of  the  normal  modes  for  the  time-harmonic  problem.  The  re¬ 
flected  froms  also  exhibit  this  feature.  The  evanescent  tails 
penetrate  deeper  into  the  water  for  these  fronts,  which  propa.- 
gate  at  higher  angles.  In  analogy,  the  length  of  the  evanescent 
tails  of  the  normal  modes  increase  with  propagation  angle.  In 
the  upslope  region,  the  evanescent  tails  grow  extremely  long 
and  energy  penetrates  into  the  sediment.  This  behavior  is  anal¬ 
ogous  to  the  coupling  of  energy  from  the  discrete  spectrum 
into  the  continuous  spectrum  for  the  time-harmonic  problem 
as  described  in  [5].  The  tails  begin  to  subside  in  length  after 
the  top  of  the  slope  is  reached. 

A  time-harmonic  source  with  u  =  100  its"'  is  approxi¬ 
mated  by 


sin  («0  =1.267  ^ 

n 


(43) 


The  constants  a"  =  1.267  •  (-1)"  and  t”  =  (n  +  jlir/w 
were  determined  by  taking  q  =  1 .  It  is  evident  from  Table  I 
that  (43)  gives  an  accurate  approximation.  By  superposition, 
the  time-harmonic  response  pcw  is  approximated  by 


Pcw(f)=  1.267  2  (-1)> 

n 


(44) 
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TABLE  I 

APPROXIMATION  OF  A  SINUSOID  BY  GAUSSIANS. 
TWO  GAUSSIANS  USED  PER  CYCLE  CENTERED 
AT  THE  MAXIMA  AND  MINIMA  OF 
THE  SINUSOID 


f(ms) 

sin  (ul) 

Xa'C.II  -  I") 

0.5 

0.1564 

0.1564 

1.0 

0.3090 

0.3089 

1.5 

0.4540 

0.4539 

2.0 

0.5878 

0.5877 

2.5 

0.7071 

0.7071 

3.0 

0.8090 

0.8091 

3.5 

0.8910 

0.8912 

4.0 

0.9511 

0.9513 

4.5 

0.9877 

0.9880 

5.0 

1.0000 

1.0003 

Fig.  2.  Benchmark  calculation  for  50  Hz  source.  Data  for  solid  curve  is 
generated  from  a  PWE  (time-domain)  calculation  followed  by  a  convolution 
to  approximate  the  CW  response.  Data  for  dashed  curve  is  generated  with 
a  PE  (frequency-domain)  calculation. 

where  p  is  the  response  to  the  Gaussian.  Fig.  2  contains  a 
plot  of  transmission  loss  at  the  receiver  depth  Zr  =  25  m 
calculated  using  (18)  with  (44)  and  using  (9).  The  excellent 
agreement  demonstrates  the  accuracy  of  the  PWE,  the  starting 
field,  and  the  numerical  solution. 

VII.  Conclusions 

A  model  that  handles  low-frequency,  bottom-interacting 
pulse  propagation  in  the  ocean  has  been  derived,  and  a  new 
initial  field  has  been  considered.  The  PWE  accounts  for  varia¬ 
tions  in  sound  speed,  density,  and  attenuation.  The  attenuation 
is  assumed  to  depend  linearly  on  frequency  in  the  sediment,  al¬ 
though  the  approach  can  be  generalized  to  handle  an  arbitraiy 
dependence.  Cylindrical  spreading  is  handled  analytically,  an 
additional  improvement  in  the  model.  A  numerical  solution 
for  the  PWE  was  derived,  and  the  accuracy  of  the  asymp¬ 
totics,  numerics,  and  starting  field  was  demonstrated  with  a 
benchmark.  The  numerical  solution  is  based  on  the  approach 
of  [17].  However,  the  solution  has  been  simplified,  and  meth¬ 
ods  designed  for  linear  problems  have  been  implemented.  The 
results  presented  here  should  apply  to  the  TDPE. 

Rirther  development  of  the  PWE  is  needed.  Future  studies 


to  extend  the  present  work  might  compare  the  PWE  with  the 
TDPE  and  compare  run-times  for  these  models  and  frequency- 
domain  models.  It  would  be  advantageous  to  extend  the  model 
to  handle  wide-angle  propagation,  which  is  known  to  be  im¬ 
portant  from  frequency-domain  studies.  The  causal  dispersion 
law  [27]  that  corresponds  to  the  linear  attenuation  law  should 
be  implemented  into  the  model  because  dispersion  is  known  to 
have  a  significant  effect  on  propagation  in  solids  [23],  and  ine 
acausal  linear  attenuation  law  sometimes  gives  unsatisfactory 
predictions  for  bottom-interacting  propagation  in  the  ocean 
[25]. 
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